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SUMMARY 

An effort is currently underway at NASA Lewis to develop two- and three-dimensional Navier-Stokes codes, 
c^ied Proteus for aerospace propulsion applications. The emphasis in the development of Proteus is 
development or research on numerical methods, but rather the development of the code itself. The <*je« u * 
develoo codes that are user-oriented, easily-modified, and well-documented. Well-proven, state-of-the-art soluuo 
algorithms are being used Code readability, documentation (both internal and external), and validauon are being 

a status report on the Proteus development effort. The analysis and soluuon procure 
rS^dbriefly -d the various Lures in the code are summarized. The results from some of the validation 
cases that have been run are presented for both the two- and three-dimensional codes. 

1. INTRODUCTION 

Much of the effort in applied computational fluid dynamics consists of modifying an existing program for what- 
J £££ and flow Isgimes arc of onrcm toes, «, che -cae«chcr. UnfaumaKly. nearly ^1 of 0* avadjNe 
nonproprietary programs were started as research projects with the emphasis on demonstrating the numerical algo- 
rithm ratheMhan SHf use or ease of modification. The developers usually intend to clean up and formally docu- 
ment the program, but the immediate need to extend it to new geometries and flow regimes takes precedence. 

The result is often a haphazard collection of poorly written code without any consistent structure. An extern 
sivelv modified program may not even perform as expected under certain combinations of operating options. Each 
new user must invent considerable time and effort in attempting to understand the underlying swclw ^^ pr ®' 
^Tintending to do anything more than run standard test cases with iL The user's subsequent modifications 
farther obscure the program structure and therefore make it even more difficult for others to understand. 

The Proteus two- and three-dimensional Navier-Stokes computer codes are intended to be urer-on«iied and 
easily-modifiable fj ow analysis programs, primarily for aerospace propulsion applications. Readabi lty, m u y, 
r y — ^ in the^primary objective. Every subroutine contain an e = ve 
describing the purpose, input variables, output variables, and calling sequence of the subrouune. Wtfo 
clearly-defined excepuons. the entire program is written in ANSI standard Fortran 77 to enhrtnce^port^ J- 
mnetw version of the program is maintained and periodically updated with corrections, as wel as ex ge - 

etal interest, such as turbulence models. 

The documentation is divided into three volumes. Volume 1 is the Analysis Description, and P^scnUthe^ua- 
Uons and solution procedure used in Proteus. It describes in detail the govemmg equations, the turbulence 
the linearization of the equations and boundary conditions, the time and space differencing formulas, ' 

viscosity Zeis. Volume 2 is the User's Guide, and contains into mauon neried 

to run the program. It describes the program’s general features, the input and output, the procedure g P 

initial conditions. the computer resource requirements, the diagnostic messages that m- 
trol language used to run the program, and several test cases. Volume 3 is the Programmer s Reference, ana con 
tainsdetailed information useful when modifying the program. It describes the program structure, the Fortran vari- 
ables stored in common blocks, and the details of each subprogram. 
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In this paper, the analysis and solution procedure are described briefly, and the various features in the code axe 
summarized. The results from some of the validation cases that have been run are presented for both the two- and 
three-dimensional codes. The paper concludes with a brief status report on the Proteus development effort, includ- 
ing the work currently underway and our future plans. 

2. ANALYSIS DESCRIPTION 

In this section, the governing equations, the numerical solution method, and the turbulence models are described 
briefly. For a much more detailed description, see Volume 1 of the documentation (Towne, Schwab, Benson, and 
Suiesh, 1990). 

2.1 GOVERNING EQUATIONS 

The basic governing equations are the compressible Navier-Stokes equations. In Cartesian coordinates, the 
two-dimensional planar equations can be written in strong conservation law form using vector notation as 1 

3Q 3E dF 3Ev dFv_ 

3 / + dx + ' = ' ' 
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The shear stresses and heat fluxes are given by 


1. pof brevity, in most instances this paper describes the two-dimensional Proteus code. The extension to three dimoisions is relatively 
straightforward. Differences between the two-dimensional and three-dimensional codes are nosed where relevant. 
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b lire* equations, t represents unte; x and , represent the Canesian coordinate directions; « tmd » « to veto- 
ties in the x and v directions- p, p, and T are the static density, pressure, and temperature, £7 is total 
rt ^Tmei n X. and X J the coefficient of viscosity, the second coefficient of viscosny. and be coefficient 

of theimal conductivity. 

In addition to the equations presented above, an equation of state is required to relate P^ejo ^"* nl 
variables. The equation currently built into the Proteus code is the equation of sate for thermally perfect gase . 
p = pRT, where R is the gas constant. For calorically perfect gases, this can be rewritten as 


/> = (/- 1) 


E T -^p(u 2 + v 2 ) 


(4) 


where y is the ratio of specific heats. c p /c v . Additional equations are also used to define p, A, k, and c„ in terms of 
temperature for the fluid under consideration. 

All of the equations have been nondimensionalized using appropriate normalizing wndiuons^ l^gths have 
been nondimensionalized by L r , velocities by u r . density by p r , temperauire by T r , viscosity by p„ erm 2 al 
Srt ^ and total energy by Pr u 2 , time by L,/u„ and gas omsunt and specific hc al by u, /T,. ^ 
reference Reynolds and Prandtl numbers are thus defined as Re, = p,u,L, I p, and Pr, - p,u,/k, ,. 

Because the governing equations are written in Cartesian coordinates, they are W -r!l; ^simplifies 
geometric configurations. For most applications a body-fitted coordinate system 

L Mjniicaiion of boundary conditions and the bookkeeping in the numerical method used to solve the equation^ 
^e duisuansformed from physical (x.y.r) coordinates to rectangular onhogonal computauonal 

coordinates. Equation (1) becomes 


dQ 3E 3F _ B¥ v 

+ d£ + drj d£ Btj 


(5) 


where 



E = -y (E£, + F£j, + Q£i) 
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F = j{Er, I + Fr, y +Qr,,) 


L-y(Ev{, + F„f,) 
FV = y(EvJ] x + FyTJy) 


In these equations the derivatives £,, r\ x , etc., are the metric scale coefficients for the generalized nonorthogonal 
grid transformation. J is the Jacobian of the transformation. 

22 NUMERICAL METHOD 

222 Time Differencing. The governing equations are solved by marching in time from some known set of initial 
conditions using a finite difference technique. The time differencing scheme currently used is the generalized 
method of Beam and Warming (1978). With this scheme, the time derivative term in equation (5) is written as 


dQ = AQ* 
dr At 


g. d(AQ") { 1 3Q" ' 

l + 0 2 dr l + 0 2 dr 


0 2 AQ 
l + 0 2 At 


- + 0 \ 
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( 6 ) 


where AQ* = Q* +1 -Q". The superscripts n and n + 1 denote the known and unknown lime levels, respectively. By 
choosing appropriate values for 0, and 0 2 , the solution procedure can be either first- or second-order accurate in 
time. 

Solving equation (5) for dQ / dr, substituting the result into equation (6) for d(AQ ") /dr and dQ /dr, and multi- 
plying by At yields 
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* n a 

2 Linearization Procedure. Equation (7) is nonlinear, since, for example, AE = E - E and the unknown 
E* is a nonlinear function of the dependent variables and of the metric coefficients resulting from the generalized 
grid transformation. The equations must therefore be linearized to be solved by the finite difference procedure. For 
the inviscid terms, and for the non-cross-derivative viscous terms, this is done by expanding each nonlinear expres- 
sion in a Taylor series in time about the known lime level n. The cross-derivative viscous terms are simply lagged 
(i.e., evaluated at the known time level n and treated as source terms.) 

The linearized form of equation (7) may be written as 
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where 3E/30 and 9F /3Q arc ihe Jacobian coefficient matrices resulting from the linearization of ** 

/SQ « *• J-*- coefficient mmriccs mridring *- *e I—— - •» 

viscous terms. 

are therefore linearized using the same procedure as for the governing equations. 

,,, solution Procedure. The governing equations. presented in linearized matrix form as equMion ' (8). are 
'2* Sgl^S^tion. equaUon (8) can be 

split into the following two-sweep sequence. 
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Sweep 2 (rj direction) 
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m* equatio ns mprcaml the two-sweep aftemnUng direcrion implicit (ADI) algorithm used to advance the solution 
from time level n to n + 1 . Q is the intermediate solution. 

Spatial derivatives in equations (9a) and (9b) are approximated using second-order central 
Hie resulting set of algebraic equations can be written in matrix form with a block tn-diagonal 1 JJ““j 

using die bksk matrix version of the THomas tdgoridun (e.g.. see Anderson. Tannehtll. and 

Pletcher. 1984). 

2.2 4 Artificial Viscosity. With the numerical algorithm described above, high frequency nonlinear instabUmes 
as die S develops. For example, in high Reynolds number flows osciUaUons can result fromthe 
odd^eruiecoupling inherent in the use of second-order central differencing for the invtscid differcnc€ 

physical phenomena such as shock waves can cause instabiliues when they are captured b * ess high 
algorithm Artificial viscosity, or smoothing, is normally added to the soluuon algorithm to suppress these high 
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frequency instabilities. Two artificial viscosity models are currently available in ti*- Proteus computer code — a 
constant coefficient model used by Sieger (1978), and the nonlinear coefficient model of Jameson, Schmidt, and 
Turkel (1981). The implementation of these models in generalized nonorthogonal coordinates is described by Pul- 
liam (1986). 

The constant coefficient model uses a combination of explicit and implicit artificial viscosity. The standard 
explicit smoothing uses fourth-order differences, and damps the high frequency nonlinear instabilities. Second- 
order explicit smoothing, while not used by Sieger or Pulliam, is also available in Proteus. It provides more 
smoothing than the fourth-order smoothing, but introduces a larger error, and is therefore not used as often. The 
implicit smoothing is second order and is intended to extend the linear stability bound of the fourth-order explicit 
smoothing. 

The explicit artificial viscosity is implemented in the numerical algorithm by adding the following terms to the 
right hand side of equation (9a) (i.e., the source term for the first ADI sweep.) 

f ^(V t A c Q + V,A,Q)--^[(V e A { ) 2 Q + (V,A,) l Q] 

ejP and cj? 5 are the second- and fourth-order explicit artificial viscosity coefficients. The symbols V and A are 
backward and forward first difference operators. 

The implicit artificial viscosity is implemented by adding the following terms to the left hand side of the equa- 
tions specified. 


-^[W'AQ)] 

--^[v,A,(7AQ")] 


to equation (9a) 
to equation (9b) 


The nonlinear coefficient artificial viscosity model is strictly explicit. Using the model as described by Pulliam 
(1986), but in the current notation, the following terms are added to the right hand side of equation (9a). 
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The subscripts i and j denote grid indices in the £ and rj directions. In the above expression, w is defined as 

¥ = ¥, + ¥, 


where Wi and w y are spectral radii defined by 


¥* = 






ivi -foV^-t-ny 

Arj 


Here U and V are the contravariant velocities without metric normalization, defined by 


and a = Vy/?7", the speed of sound. 


U - + 4 iU + {yV 

V = r], + rhu + T] y v 
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■tte r .... . £ ® and t" are the second- and fonnh-order «ii 6 cial .iscosil, coefficiems. For d» coefecienis 

of the 4 direction differences, 


JejPj = K- 2 A*max(a, +1 , a,, <r,.\) 
^ 4) j = max 0 , k a At- [ £ ? > ] . 


where 




\ pM-1Pi+Pi-l | 

|Pi + l + 2/>i + Pi-l | 


and r, and *4 are constants. Similar formulas are used for the coefficients of the rj direct 10 " TJ* 

^Lneter o is a pressure gradient scaling parameter that increases the amount of 

tofourth -order smoothing near shock waves. The logic used to compute e< switches off the fourth-order smooth 
ing when the second-order smoothing term is large. 


13 TURBULENCE MODELS 

Turbulence is modeled using either a generalized version of the Baldwin and Lomax (1978) algebraic eddy 
viscosity model, or the Chien (1982) low Reynolds number k-c model. 

2J.1 Baldwin-Lomax Model. For wall-bounded flows, the Baldwin-Lomax turbulence model is a two-layer 
model, with 


fO*.W toy K <y b (10 ) 

|OOou*r for y K >y b 


where y, is the normal distance from the wall, and y b is the smallest value ofy. at which the values of p, from the 
inner and outer region formulas are equal. For free turbulent flows, only the outer region value is used. 


The outer region turbulent viscosity at a given 4 or r? station is computed from 

( Pi)outtr = KC cpPF KUb^ wakt^r 


(ID 


where ms the Clauser constant, taken as 0.0168, and C cp is a constant taken as 1.6. 
The parameter is computed from 


ymax^max for wall-bounded flows 

for free turbulent flows 

F mm 


( 12 ) 


where C** is a constant taken as 0.25, and 

Vjiff = 1^1 max~ 1 ^ 1 «« 


where is the total velocity vector. 

The parameter in equation (12) is the maximum value of 
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F(y„) 




-y'M 
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for wall-bounded flows 
for free turbulent flows 


(13) 


and v is the value of y„ corresponding to F^. 

Fbr wall-bounded flows, y„ is the normal distance from the wall. For free turbulent flows, two values of F^ 

and y are computed — one using the location of I V as the origin for y m , and one using the location of 

,T? Uv The origin giving the smaller value of y„, is the one finally used for computing y„, F mua , and y*. 

In equation (13), |fl | is the magnitude of the total vorticity, defined for two-dimensional planar flow as 

till = < l4 > 

| OX oy | 

The parameter A * is the Van Driest damping constant, taken as 26.0. The coordinate y + is defined as 

* 1 „ 'J^wPwFCr /IS! 

y* — Re,= y„ (>31 

Pw Pw 


where u, = -fajpjie, is the friction velocity, r is the shear stress, and the subscript w indicates a wall value. In 
Proteus, t w is set equal to/r,JflL. 

The function F K ub in equation (11) is the Klebanoff intermiltency factor. Fbr free turbulent flows, F KUb = 1. 
For wall-bounded flows. 


Fxub - 



61-1 


(16) 


In equation (16), B and C Kub are constants taken as 5.5 and 0.3, respectively. 
The inner region turbulent viscosity in the Baldwin-Lomax model is 

OOi««r = p / 2 l^lfor 


(17) 


where l is the mixing length, given by 


(18) 


and r is the Von Karman constant, taken as 0.4. 

If both boundaries in a given coordinate direction are solid surfaces, the turbulence model is applied separately 
for each surface. An averaging procedure is used to combine the resulting two n, profiles into one. 

The turbulent second coefficient of viscosity is simply defined as 



The turbulent thermal conductivity coefficient is defined using Reynolds analogy as 


c p H, 

Pr 


where c p is the specific heat at constant pressure, and Pr, is the turbulent Prandil number. 


138 



kinetic energy and ihe turbulent dissipation rate, respectively. 

„ cmsian Mortis. te .wo-dimenrlonai planar aquadom for UK Chian is modal can ba whuan urn, 

vector notation as 
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The turbulent viscosity is given by 


(20) 


In the above equations, C,. C^, C 3 , <x t , <r r . and C„, are constants equal to 1.35, 1.8, 0.01 15. 1.0. 1.3. and 0.09, 
respectively. The parameter y n is the minimum distance to the nearest solid surface, and y is computed from y „ . 
In the above equations the mean flow properties have been nondimensionalized as described in Section 2.1. The 
turbulent kinetic energy k and the turbulent dissipation rate e have been nondimensionalized by w? and p,u, /p r , 
respectively. 

After transforming from physical to rectangular orthogonal computational coordinates, equation (19) becomes 
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The time differencing scheme and linearization procedure described previously for the mean flow equations are 
also applied to equation (21). The mean flow variables are evaluated at the known tune level n. This ** 
equations to be uncoupled from the mean flow equations and solved separately^ Spatid J^J^ffere^ es foMhe 
mated using first-order upwind differences for the convective terms, and second-order central ‘ Z fe 

viscous terms. In the two-dimensional Proteus code, the equations are solved by the same ADI procedure as the 
mean flow equations. In the three-dimensional code, they are solved by a two-sweep LU procedure, as described by 

Hoffmann (1989). 

The turbulent second coefficient of viscosity A, and the turbulent thermal conductivity coefficient k, are defined 
as described in the previous section. 

3. CODE FEATURES 

In this section the basic characteristics and capabilities of the two- and three-dimensional Proteus cries i are 
summarized. For a much more detailed description, see Volumes 2 and 3 of the documentation (Towne, Sch a . 
Benson, and Suresh, 1990). 

3.1 ANALYSIS 

The Proteus codes solve the unsteady compressible Navier-Stokes equations in either two or three dimensions^ 
The 2-D code can solve either the planar or axisymmetric form of the equations. Swirl is allowed m 
flow. The 2-D planar equations and the 3-D equations are solved in fully conservative form. As subsets of these 
equations, options are available to solve the Euler equations or the thin-layer Navier-Stokes equations. An op 
also available to eliminate the energy equation by assuming constant total enthalpy. 

Hie equations are solved by marching in time using the generalized time differencing of Beam and Earning 
(1978). The method may be either first- or second-order accurate in time, depending on the choice c rf ume ** 
(denting parameters. Second-order central differencing is used for all spatial denvauves. Nonlinear tenns are 
linearized Sing second-order Taylor series expansions in time. Hie resulting difference equations are solved using 
an alternating-direction implicit (ADI) technique, with Douglas-Gunn type splitting as written by Bn ey an 
McDonald (1977). The boundary conditions are also treated implicitly. 

Artificial viscosity, or smoothing, is normally added to the solution algorithm to damp pre- and post-shock oscil- 
lations in supersonic flow, and to prevent odd-even decoupling due to the use of central differences in cmvecuon- 
dominated regions of the flow. Implicit smoothing and two types of explicit smoothing are available in ■ 

Hie implicit smoothing is second order with constant coefficients. For the explicit smoothing the 
constant coefficient second- and/or fourth-order model (Sieger, 1978), or a nonlinear <^ffic.ent mixed «cand- and 
fourth -order model (Jameson, Schmidt, and Turkel, 1981). The nonlinear coefficient model was designed 
specifically for flow with shock waves. 

The equations are fully coupled, leading to a system of equations with a block tridiagonal coefficient matrix that 
can be solved using the block matrix version of the Thomas algorithm. Because this algorithm is recursive the 
source code cannot be vectorized in the ADI sweep direction. However, it is vectorized in the non-sweep direction, 
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leading to an efficient implementation of (he algorithm. 

3 2 GEOMETRY AND GRID SYSTEM 

The equations solved in Proteus were originally written in a Cartesian coordinate system, then transformed into 
a general nonorthogonal computational coordinate system. The code is therefore not limited to any particular type 
of geometry or coordinate system. The only requirement is that body-fitted coordinates must be used. In general, 
the computational coordinate system for a particular geometry must be created by a separate coordinate generation 
code and stored in an unformatted file that Proteus can read. However, simple Cartesian and polar coordinate sys- 
tems are built in. 

The equations are solved at grid points that form a computational mesh within this computational coordinate 
system. The number of grid points in each direction in the computational mesh is specified by the user. The loca- 
tion of these grid points can be varied by packing them at either or both boundaries in any coordinate direction. The 
transformation metrics and Jacobian are computed using finite differences in a manner consistent with the differenc- 
ing of the governing equations. 

3 3 FLOW AND REFERENCE CONDITIONS 

As stated earlier, the equations solved by Proteus are for compressible flow. Incompressible conditions can be 
simulated by running at a Mach number of around O.L Lower Mach numbers may lead to numerical problems. The 
flow can be laminar or turbulent. The gas constant R is specified by the user, with the value for air as the default 
The specific heats c p and c„, the molecular viscosity p, and the thermal conductivity A: can be treated as constants or 
as functions of temperature. The empirical formulas used to relate these properties to temperature ate contained in a 
separate subroutine, and can easily be modified if necessary. The perfect gas equation of state is used to relate pres- 
sure. density, and temperature. This equation is also contained in a separate subroutine, which could be easily 
modified if necessary. All equations and variables in the program are nondimensionalized by normalizing values 
derived from reference conditions specified by the user, with values for sea level air as the default. 

3.4 BOUNDARY CONDITIONS 

The easiest way to specify boundary conditions in Proteus is by specifying the type of boundary (e.g., no-slip 
adiabatic wall, subsonic inflow, periodic, etc.). The program will then select an appropriate set of conditions for that 
boundary. For many applications this method should be sufficient. If necessary, however, the user may instead set 
the individual boundary conditions on any or all of the computational boundaries. 

A variety of individual boundary conditions are built into the Proteus code, including: (1) specified values 
and/or gradients of Cartesian velocities u, v, and w, normal and tangential velocities V„ and V„ pressure p, tempera- 
ture T, and density p\ (2) specified values of total pressure pr, total temperature T t , and flow angle; and (3) linear 
extrapolation Another useful boundary condition is a "no change from initial condition" option for u, v, w, p, T, p, 
Pr, and/or T t . Provision is also made for user-written boundary conditions. Specified gradient boundary conditions 
may be in the direction of the coordinate line imersecting the boundary or normal to the boundary, and may be com- 
puted using two-point or three-point difference formulas. For all of these conditions, the same type and value may 
be applied over the entire boundary surface, or a point-by-point distribution may be specified. Unsteady and time- 
periodic boundary conditions are allowed when applied over the entire boundary. 

3.5 INITIAL CONDITIONS 

Initial conditions are required throughout the flow field to start the time matching procedure. For unsteady flows 
they should represent a real flow field. A converged steady-state solution from a previous run would be a good 
choice. For steady flows, the ideal initial conditions would represent a real flow field that is close to the expected 
final solution. 

The best choice for initial conditions, therefore, will vary from problem to problem. For this reason Proteus 
does not include a general-purpose routine for setting up initial conditions. The user must supply a subroutine, 
called INIT, that sets up the initial starting conditions for the time marching procedure. A version of INTT is, how- 
ever, built into Proteus that specifies uniform flow with constant flow properties everywhere in the flow field. These 
conditions, of course, do represent a solution to the governing equations, and for many problems may help minimize 
starting transients in the time marching procedure. However, realistic initial conditions that are closer to the 
expected final solution should lead to quicker convergence. 
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3.6 TIME STEP SELECTION 

Semal differed opuons ait available for choosing 0* lime amp 4.. and to moJilyi^i; as *u soUmoo 

sz!s2xsgsggs£35&£&£ 

rSiSSss^r ; 

available to cycle At between two values in a logarithmic progression over a specified number of ume P 

3.7 CONVERGENCE 

Five options are currently available for determining convergence. The user specifies a cOT^g^ c^nOTf 

It should be noted however, that convergence is in the eye of the beholder. The amount of decrease ^ 
dual necessary for convergence will vary ftom problem to problem. For some problems, it may ^ m ° re **£ ? . 
il^Snv^nce by some floriated parameter, such as the lift coeffimtt for an airfo.1. Determ, rung 
whenT^lution is sufficiently converged is, in some respects, a skill best acquired through experience. 

3.8 INPUT/OUTPUT 

Incut to Proteus is through a series of namelists and, in general, an unformatted file containing the computa- 

tion/coordinaie system. All of the input parameters have default values ‘ l ° * ^ “ “n Is 

a different value is desired. Reference conditions may be specified in either English or SI units. A P 

also available, in which the computational mesh and the initial flow field are read from un ortna 

created during an earlier run. 

The standard printed output available in Proteus includes an echo of the input, boundary conditions, normalizuig 
computed flow field, end conveigence infomudion. TV use- condole 
kL field parameters are printed, and at which time levels and grid points. Several debug opuons are also avadable 
for detailed printout in various parts of the program. 

In addition to the printed output, several unfotmatted files can be written for various pwposes. ****** 
auxiliary file used for post-processing, usually called a plot file, that can be written at convergence or after the last 
time sum if the solution docs not converge. Plot files can be written for the NASA Lewis piotung progr^ CON- 
TOURot the NASA Ames plotting program PLOt3D. If PLOT3D is to be used, two unformatted files are created, 
an xyz file containing the computational mesh and a q file containing the computed flow field Anoiher unformat 

Proteus contains detailed convergence information. This file is automatically incremented each i time 
Se Sn is checked for convergence, and is used to generate the convince 

developed post processing plotting routines. And finally, two unformatted files may be written at the end of a 
lationthat may be used to restart the calculation in a later run. One of these contains the computational mesh, an 
the other the computed flow field. 

3.9 TURBULENCE MODELS 

For turbulent How. Fromm -ol.es Ihu Reynolds Umo-uverngod N«.ie,-S t ok« e,m^ns_ wUh uubulence 
modeled using eilher die Baldwin and Umax (1978) algebraic eddyviscosuy model or die Chien (1982) two- 
equation model. 

33.1 Baldwin-Lomax Model. The Baldwin-Lomax model may be applied to either wall-bounded flows or tofree 
turbulent flows. For wall-bounded flows, the model is a two-layer model. For flows in which more than one boun 
dary is a solid surface, averaging procedures are used to determine a single n, profile. The turbulent thermal 
ductivity coefficient k, is computed using Reynolds analogy. 

3 9 2 Chien j k-e Model. With the Chien two-equation model, partial differential equations are solved for the tur- 
bint kmetic enerey * and the turbulent dissipation rate «. These equations are lagged in time and solved 
5 "ow equations. In the 2-D Proteus code, the equations are solved using the same solution 
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algorithm as for the mean flow equations, except that spatial derivatives for the convective terms are approximated 
using first-order upwind differencing. In the 3-D code, they are solved by a two-sweep LU procedure, as described 
by Hoffmann (1989). 

Since the Chien two-equation model is a low Reynolds number formulation, the k-e equations are solved in the 
near-wall region. No additional approximations are needed. Boundary conditions that may be used include: (1) no 
r hang* from initial or restart conditions for k and e\ (2) specified values and/or gradients of k and e; and (3) linear 
extrapolation. Specified gradient boundary conditions are in the direction of the coordinate line intersecting the 
boundary, and may be computed using two-point or three-point difference formulas. For all of these conditions, the 
same type and value may be applied over the entire boundary surface, or a point-by-point distribution may be 
specified. Spatially periodic boundary conditions for k and e may also be used. Unsteady boundary conditions are 
not available for the k-e equations. However, unsteady flows can still be computed with the Chien model using the 
unsteady boundary condition option for the mean flow quantities and appropriate boundary conditions for k and e, 
such as specified gradients or linear extrapolation. 

initial conditions for k and e are required throughout the flow field to start the time marching procedure. The 
best choice for initial conditions will vary from problem to problem, and the user may supply a subroutine, called 
KEINIT, that sets up the initial values of k and e for the time marching procedure. A version of KEINIT is built 
into Proteus that computes the initial values from a mean initial or restart flow field based on the assumption of local 
equilibrium (i.e., production equals dissipation.) Variations of that scheme have been found to be useful in comput- 
ing initial k and e values for a variety of turbulent flows. 

The time step used in the solution of the k-e equations is normally the same as the time step used for the mean 
flow equations. However, the user can alter the time step, making it larger or smaller than the time step for the 
mean flow equations, by specifying a multiplication factor. The user can also specify the number of k-e iterations 
per mean flow iteration. 


4. VERIFICATION CASES 

Throughout the Proteus development effort, verification of the code has been emphasized. A variety of cases 
have been run, and the computed results have been compared with both experimental data and exact solutions. 
Some cases are included in Volume 2 of the Proteus documentation (Towne, Schwab, Benson, and Suresh, 1990). 
Other cases have been reported by Conley and Zeman (1991), Saunders and Keith (1991), and Bui (1992). 

Three cases are presented in this paper — flow past a circular cylinder, flow through a transonic diffuser, and 
flow through a square-cross-sectioned S-duct. 

4.1 FLOW PAST A CIRCULAR CYLINDER 

In this test case, steady flow past a two-dimensional circular cylinder was investigated. Both Euler and laminar 
viscous flow were computed. 

4.1.1 Reference Conditions. In order to allow comparison of the Proteus results with incompressible experimental 
and with potential flow results, this case was run with a low reference Mach number of 0.2. The cylinder 

radius was used as the reference length, and was set equal to 1 ft Standard sea level conditions of S19 °R and 
0.07645 lbm / ft 3 were used for the reference temperature and density. The Reynolds number based on cylinder 
diameter was 40, matching the experimental value. 

4.12 Computational Coordinates. For this problem a polar computational coordinate system was the obvious 
choice. The radial coordinate r varied from 1 at the cylinder surface to 30 at the outer boundary. Since the flow is 
symmetric, only the top half of the flow field was computed. The circumferential coordinate 6 thus varied from 0° 
at the cylinder leading edge to 180° at the trailing edge. For the Euler flow case, a 21 (circumferential) x 51 (radia 1 ) 
mesh was used, with the radial grid packed moderately tightly near the cylinder surface. For the viscous flow case, 
a 51 x51 mesh was used, with the radial grid packed more tightly near the cylinder surface. 

4.1.3 Initial Conditions. Constant stagnation enthalpy was assumed, so only three initial conditions were required. 
For the Euler flow case, uniform flow with u = 1, v = 0, and p = 1 was used. 

For the viscous flow case, the exact potential flow solution was used to set the initial conditions at all the non- 
wall points. Thus, with nondimensional free stream conditions of p m = u„ = T» = p- = 1, the initial conditions 
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were 2 

U = 1 rCOS(20) 

r 


v = - -YSin(20) 

r 

1 p„(u 2 + v 2 ) 
P = (Pt)-~ 2 p 


where 

l 

(p T )- =P-+2 

surface 

4.1.4 Boundary Conditions. Again, since. 

lions were required ai each computational boundap'. «.rface the radial velocity and the radial gra- 

the symmetry line ahead of and behind the cy in er i t y di ^ l of pressure was computed from the 

diem of the circumferential velocity were set equal to zero. The rndal oM»e s v 

polar coordinate form of the incompressible rad.al momentum equauon wnlten at the watt, t ne eq 

and Gaylord, 1964) 

dv, pv e dv r v[ dp 

P V '-dT+-rjr P r "dr 

where v r and v, are the radial and circumferential velocities, respectively. At the cylinder surface, v, = 0. Thus, 


d/L = 

dr 


Ve 

P—=P 


u 2 +v 2 


And finally, at the outer boundary the tree stream conditions were specified as boundary rood, uons. 

For ihe viscous flow case, symmeiry ctmdiiions were gradient 

S"^ren!i” flol" ^InT".^ SSeUn^ ^loes of , were hep. a. dreir initial 
values, and the radial gradients of u and v were set equal to zero. 

4>1J Numerics. Both the Euler and viscous flow cases were run t*ing a ^ 

CFL number of 10. The constant coefficient artificial viscosity model was used, wi , i 

The Euler flow case converged in 2l0ume asps , and ^ ^ 

convergence criterion for both cases was that the L 2 norm of * c rcsiclual OT ^ 


tv.. V*,..,,,. in the Proteus inpul ind output, the pressure is 

2. Note lh»t the nondunension.l g*s consum R »ppe*n in these iqMM described m Section 2 1 
nondimensiontlized by pJtT,. Intemtl to the code, pressure is nondurtenstonsb^ by p.u , . 
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4.1.6 Computed Results. In Figure 1 the computed static pressure coefficient, defined as (p-p r )/(Pr^ /2g e ) is 
plotted as a function of 6 for both the Euler and viscous flow cases. Also shown are the experimental data of Grove, 
Shair. Petersen, and Acrivos (1964), and the exact solution for potential flow. The Proteus results agree well with 
the data for the viscous flow case, and with the exact potential flow solution for the Euler flow case. 



Figure 1. Pressure coefficient for flow past a circular cylinder. 


4.2 TRANSONIC DIFFUSER FLOW 

In this test case, two-dimensional transonic turbulent flow was computed in a converging-diverging duct. Tur- 
bulence was modeled using the Baldwin-Lomax model. The flow entered the duct subsonically, accelerated through 
the throat to supersonic speed, then decelerated through a normal shock and exited the duct subsonically. The com- 
putational domain is shown in Figure 2. 


±- 



Figure 2. Computational domain for transonic diffuser flow. 


4.2.1 Reference Conditions. The throat height of 0.14435 ft. was used as the reference length L r . The reference 
velocity u, was 100 ft/sec. The reference temperature and density were 525.602 °R and 0.1005 lb,* /ft , respec- 
tively. These values match the inlet total temperature and total pressure used in other numerical simulations of this 
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flow (Hsich, Bogar, and Coakley, 1987). 

4 2 2 Computational Coordinates. The x coordinate for this duct runs from -4.04 io +8- 6 5. The Cartes^ coordi _ 
nates of the bottom wall are simply y = 0 for all x. For the top wall, the y coordinate is given by (Bogar. Sajben, and 

Kroulil. 1983) 


1.4144 

y =^acoshf /(a-l+coshf) 
1.5 


for x £ -2.598 

for -2.598 <x <7.216 

for x £7.216 


where the parameter f is defined as 


Ci(x/x t )[l + C 2 x/x,] C ’ 
(1 -x/x,) C * 


The various constants used in the formula for the top wall height in the converging (-2.598 <x SO) and diverging 
(OSx S7.216) parts of the duct are given in the following table. 


Constant 

Converging 

Diverging 

a 

1.4114 

1.5 

x, 

-2.598 

7.216 

c, 

0.81 

2.25 

c 2 

1.0 

0.0 

c 3 

0.5 

0.0 

c 4 

0.6 

0.0 


A body-fitted coordinate system was generated for the duct, with 81 points in the x direction and 51 pontu “t «he 
y direction The coordinate system is shown in Figure 2. For clarity, the grid points are thinned by facwrs of 2 and 
10 in the x and y directions, respectively. Note that for good resolution of the flow near the normal shock, the gnd 
defining the computational coordinate system is denser in the x direction m the region just downst^ oft^thrM 
In the y direction, the actual computational mesh was tightly packed near both walls to resolve the turbulent boun- 
dary layers. 

4 2. 3 Initial Conditions. The initial conditions were simply zero velocity and constant pressure and temperature. 
Thus, u = v = 0 and p = T= 1 everywhere in the flow field. 

4.2.4 Boundary Conditions. This calculation was performed in three separate runs. In the first run. the exit stauc 
pressure was gradually lowered to a value low enough to establish supersonic flow throughout the diverging portion 
of the duct. The pressure was lowered as follows: 


p(0 = 


0.99 

- 2 . 1405xl(T 3 n + 1.20405 
0.1338 


for 1 <n£l00 
for 101 £n £500 
for 5015a £3001 


where a is the time level. The equation forp for 101 £n £500 is simply a linear interpolation between p = 0.99 and 
p = 0.1338. In the second run. the exit pressure was gradually raised to a value consistent with the formauon o a 
normal shock just downstream of the throat Thus, 

f3.4327xl0^n-0.89636 for 3001 < n £ 5000 


p(0 = 


0.82 


for 5001 <n 56001 


Again, the equation for p for 3001 < n < 5000 is simply a linear interpolation between p - 0.1338 and p - 0.82. In 
the third run, the exit pressure was kept constant at 0.82 for 6001 < n < 9000. 
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The remaining boundary conditions were the same for all runs. At the inlet, the total pressure and total tempera- 
ture were set equal to 1, and the y-velocity and the normal gradient of the x - velocity were both set equal to zero. At 
the exit, the normal gradients of temperature and both velocity components were set equal to zero. At both walls, 
no-slip adiabatic conditions were used, and the normal pressure gradient was set equal to zero. 

4 7 5 Numerics. The case was run using a spatially varying time step. The local CFL number was 0.5 for the first 
two runs, and 5.0 for the third run. The nonlinear coefficient artificial viscosity model was used. For the first two 
runs, the coefficients c (2) and were 0.1 and 0.005, respectively. For the third run, c (4) was lowered to 0.0004. 

The convergence criterion was that the absolute value of the maximum change in the conservation variables 
AQmb be less than 10 -6 . At the end of the third run, die solution had not yet converged to this level. However, 
close examination of several parameters near the end of the calculation indicates that the solution is no longer 
changing appreciably with time, but oscillates slightly about some mean steady level. This type of result appears id 
be fairly common, especially for flows with shock waves. The reason is not entirely clear, but may be related ® 
inadequate mesh resolution, discontinuities in metric information, etc. For this particular case, the cause may also 
be inherent unsteadiness in the flow. The experimental data for this duct show a self-sustained oscillation of the 
normal shock at Mach numbers greater than about 1.3 (Bogar, Sajben, and Kroulil, 1983). 

4 . 2.6 Computed Results. The computed flow field is shown in Figure 3 in the form of constant Mach number con- 
tours. 



Figure 3. Computed Mach number contours for transonic diffuser flow. 

The flow enters the duct at about M = 0.46, accelerates to just under M = 1.3 slighdy downstream of the throat, 
shnrlrQ down to about M - 0.78, then decelerates and leaves the duct at about M = 0.51. The normal shock in the 
throat region and the growing boundary layers in the diverging section can be seen clearly. Because this is a shock 
capturing analysis, the normal shock is smeared in the streamwise direction. 

The computed distribution of the static pressure ratio along the top and bottom walls is compared with experi- 
mental data (Hsieh, Wardlaw, Collins, and Coakley, 1987) in Figure 4. The static pressure ratio is here defined as 
P /(Pr) o, where (pj-)o is die inlet core total pressure. 
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The computed results generally agree well with the experimental 
normal shock. The predicted shock position however, is s lg tyj distancc There is also some disagreement 

m ^ - 

contour in this region without sufficient mesh resolution. 

43 TURBULENT S-DUCT FLOW 

In this lea case, three-dimensional turbulent flow in an Mu' 22.5 s bends wilh 

si,om . - — * 

Taylor. Whitelaw, and Yianneskis ( 1982 ). 

4J.1 Reference Conditions. The default standajd sea level ^iuonsforjur 

used for the reference temperature and density. The specific beat ratio IT compressibility effects and. at the 

incompressible, the refetence Mach number M, »t set tepid to °f “TT” ‘SL tte Reynolds number 
same time, achieve a reasonable convergence rate «th tte Promos code. !*ed as the reference 

based on die bulk velocity and die hydraulic diame r w , - ^ q 028658 ft. This value was com- 

Reynolds number Re r in the calculation. The reference length L r was respectively 

25 from the definition of Re„ where Af r and Sutherland’s law were used to compute u r and ^ r . respecuve.y . 

4.U Computational Coordinates. Figum 5 _ . t 

GR1DCEN codes (Steinbrenner. Chawner, and (actor rftwTtn each direction. The boundary 

^^r^nTdfoSrTO 6 ^^ U ~ - - i— - - 

boundary grids using GRIDGEN 3D. 
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Figure 5. S-duct computational grid. 


The computational grid extended from 7.5 hydraulic diameters upstream of the start of the first bend, to 7.5 
hydraulic diameters downstream of the end of the second bend. The grid consisted of 81 x 31 x61 points in the £, r/, 
and f directions, respectively. Since the S-duct is symmetric with respect to the r} = 1 plane, only half of the duct 
was computed. To resolve the viscous layers, grid points were tightly packed near the solid walls using the default 
packing option in GRIDGEN 2D. At the grid point nearest the wall, the value of y * was about 0.5. 

433 Initial Conditions. The computations were done in two separate major steps: a calculation using the 
Baldwin-Lomax turbulence model and a calculation using the Chien k-c model. To start the Baldwin-Lomax calcu- 
lations, the default initial profiles specified in subroutine INIT were used. Thus, the static pressure p was set equal 
to 1.0, and the velocity components u, v, and w were set equal to 0.0 everywhere in the duct. To start the Chien k-e 
calculations, the initial values of u> v, w, p, and the turbulent viscosity p t were obtained from the Baldwin-Lomax 
solution. The initial values of k and £ were obtained using the default KEINIT subroutine in Proteus . 

43-4 Boundary Conditions. For both calculations, constant stagnation enthalpy was assumed, eliminating the 
need for solving the energy equation. Therefore, only four boundary conditions were required for the mean flow at 
each computational boundary. In addition, for the Chien calculation, boundary conditions were required for k and e 
at each computational boundary. 

For the Baldwin-Lomax calculation, at the duct inlet the total pressure was specified as 1.02828, the gradient of 
u was set equal to zero, and the velocities v and w were set equal to zero. The inlet total pressure was calculated 
from the freestream static pressure and the reference Mach number using isentropic relations. At the duct exit, the 
static pressure was specified as 0.98416, and the gradients of u, v, and w were set equal to zero. The exit static pres- 
sure was found by trial and error in order to match the experimental mass flow rate. At the walls of the duct no-slip 
conditions were used for the velocities, and the normal pressure gradient was set to zero. Symmetry conditions 
were used in the symmetry plane. 

For the Chien calculation, the boundary conditions for the mean flow were the same as for the Baldwin-Lomax 
calculation, with one exception. At the duct exit, the value of the static pressure was changed slightly, from 0.98416 
to 0.98474, again in order to match the experimental mass flow rate. For the J fc-c equations, at the upstream boun- 
dary the gradients of the turbulent kinetic energy k and the turbulent dissipation rate £ were set equal to zero for the 
first 20 time steps. After that time, the values of k and e were kept constant At the downstream boundary, the gra- 
dients of k and c were set equal to zero. No-slip conditions were used at the solid boundaries, and symmetry condi- 
tions were used at the symmetry boundary. 

433 Numerics. Both the Baldwin-Lomax and Chien calculations were run using a spatially varying time step. 
Since the flow field for the Baldwin-Lomax calculation was impulsively started from zero velocity everywhere, 
large CFL numbers specified at the very beginning of the calculation might result in an unphysical flow field and 
cause the calculation to blow up. Therefore, the calculations were run with a CFL number of 1 for the first 100 
iterations, 5 for the next 200 iterations, and 10 for the remaining iterations. A total of 4,000 iterations was used for 
the Baldwin-Lomax calculation. 
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iterations. A total of 2.520 iterations was used for the Chien calculation. 

The constant coefficient artificial viscosity model was used for both cases, with e, - 2 and 4 - 1. 

Thu criterion was that the aventge residoa! for each 

ctdations were stopped before reaching Otis level of convergence - » 5e e^7*e 

ters indicated that the solution was no longer changing appreciably «h » ' i xl( r?for the continuity equation. 
Baldwin-Lomax calculation ranged from 10* for the x-momentum equauon “ ,2 Z ^“Ltiyequa- 

Por the Chien calculation the values were 3x10- for the x-momentum equauon rtMT for the conunuity equa 
lion For both cases the residuals were continuing to drop when the calculauons were stopped. 

a i a rnmnuted Results In Figure 6. the computed flow field from the Chien calculation is shown in the fonn of 

S?flow vortices in the upper half of the cross-section. TTiese secondary flows cause a stgmficant amount of flow 
distortion, as shown by the total pressure contours. 

In the second bend, the direction of the cross-flow pressure gradients reverses, making the pressure Wghcr inthe 
„pp£ However, ihe flow emers ihe recond bend wirh . vonex pe rrem Hre«l y eri^hshed. 

The net effect is to tighten and concentrate the existing vortices near the top of the duct, in a 8 rcc ”\ . . . r 

£^“ow ibe^y. The resulting horseshoe-shap* distonion pauen. a. dre exiiof.be second bend ,s typical of 

S-duct flows. 



Figure 6. Computed total pressure contours for turbulent S-duct flow. 

In Figure 7 the calculated wall pressure distribution is compared with the experimental data of Taylor, Whi- 
relaw.^"kis(1982). The agreement is very good. Both turbulence .«*b 
sure trend and the pressure loss along the duct. The r and : coordinates noted in the legend are the same as 
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defined by Taylor, Whiielaw, and Yianneskis. 



Streamwise Distance, D h 


Figure 7. Computed surface static pressure distribution for turbulent S-duct flow. 


In Figure 8, the experimental and computed velocity profiles in the symmetry plane are shown for the five 
streamwise stations that were surveyed in the experiment. These survey stations are at the same locations as the 
total pressure contours shown in Figure 6. The agreement between computation and experiment is excellent for 
both turbulence models. The asymmetry in the velocity profiles due to the pressure induced secondary motion is 
correctly predicted. 
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Figure 8. Computed stream wise velocity profiles for turbulent S-duct flow. 


5. CONCLUDING REMARKS 

-tot Proteus two- and three-dimensional Navier-Stokes c^^ n ^vel^at 

described, and results have been presented r0 ™ ^ mc ° . Suresh 1990) and version 2.0 was released in late 

code released in tare 1989 (Towne. Schwab. Benson for version 2.0 of Ore 

- — *-» 

developmem work on dre Proreus codes is being done ioadd a nruHrple -rone grid apabiliry. a muhi- 
grid convergence acceleration capability, and additional turbulence modeling options. 

A wide variety of validation cases have been run, including: (1) sev«nl 

Navier-Stokes solutions exist; (2) laminar and ^ k wavCS . (5) steady and unsteady flows past 

dimensional driven cavity flows. (4) flows with normal and ducU5; 

a cylinder. (6) developing laminar and turbulent flows J n J*^ ^ d (9) ^rbulent flow on a flat plate 
~ C f r ; duct flows and 

flows with heat transfer. 
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